﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using ODESolver;

namespace Crystal
{
    public delegate StatusWindow stWindow();

    public class Atom
    {
        public const string ATOM_SORT_NA  =  "Na";
        public const string ATOM_SORT_I   =  "I";
        public const string ATOM_SORT_H   =  "H";

        public int deg;
        public string sort;
        public double mass;
        public double charge;
        public double[] pos;
        public double[] vel;
        /// <summary>
        /// It should be object references to atoms
        /// </summary>
        public List<object> AtomsToInteractWith;

        public InteractionDoubleWell potE 
        {
            set { _potE = value; }
            get
            {
                if (_potE == null)
                    return (dist, pars) => { return 0.0; };
                else
                    return _potE;
            }
        }
        private InteractionDoubleWell _potE;

        public string msg
        {
            set{ _stWindow().addText(value, StatusWindow.TXT_WARNING);}
        }

        /// <summary>
        /// The form that will recieve status information and bring it to user
        /// </summary>        
        private static stWindow _stWindow;
        public static void init(stWindow stWin)
        {
            _stWindow = stWin;
        }

        public Atom(double[] pos, double[] vel, string sort)
        {
            if (pos.Length != vel.Length) msg = "The vel have different size with pos";
            this.deg = pos.Length;
            this.pos = pos;
            this.vel = vel;
            this.sort = sort;
            switch (sort)
            {
                case "Na": mass = 23; charge = 11; break;
                case "I": mass = 128; charge = 53; break;
                case "H": mass = 1; charge = 1; break;
                default: mass = 1; charge = 1; break;
            }
            AtomsToInteractWith = new List<object>();
        }

        public double KineticEnergy()
        {
            double KinEnergy = 0.0;
            for (int i = 0; i < deg; i++)
            {
                KinEnergy += vel[i] * vel[i];
            }
            KinEnergy *= mass / 2;
            return KinEnergy;
        }

        public void MoveTo(double[] vect)
        {
            for (int i = 0; i < deg; i++)
            {
                this.pos[i] += vect[i];
            }

        }


        public double distance(Atom a)
        {
            if (deg != a.deg)
            {
                msg = "Can't culculate the distance of two atoms with different degries";
                return 0.0;
            }
            double dist = 0.0;
            for (int i = 0; i < deg; i++)
                dist += (this.pos[i] - a.pos[i]) * (this.pos[i] - a.pos[i]);
            return Math.Sqrt(dist);
        }

        public double Distance(double[] pos)
        {
            if (pos.Length != deg) msg = "The pos1 have different size with pos2";
            double dist = 0.0;
            for (int i = 0; i < deg; i++)
            {
                dist += (this.pos[i] - pos[i]) * (this.pos[i] - pos[i]);
            }
            return Math.Sqrt(dist);
        }

        public double Potential(double dist, string sort)
        {
            double[] potPars = SystemEquetions.getMorzePars(this.sort, sort);
            //the potencial function
            return potE(dist, potPars);
        }

        public double PairInterection(double[] pos, string sort)
        {
            double pot = Potential(Distance(pos), sort);
            pot = Math.Max(pot, -200);
            pot = Math.Min(0.5, pot);
            return pot;
        }

        public double PairInterection(Atom a)
        {
            double pot = Potential(this.distance(a), a.sort);
            pot = Math.Max(pot, -200);
            pot = Math.Min(0.5, pot);
            return pot;
        }
    }

    public class MechSys : List<Atom>
    {
        /// <summary>
        /// The form that will recieve status information and bring it to user
        /// </summary>        
        private static stWindow _stWindow;
        private SystemEquetions _eq;
        private sysTypes _interType;

        private double _currT = 0;
        private double _dt = 0.1;

        public sysTypes interType 
        {
            get { return _interType; }
            set 
            {
                _interType = value;
                _eq = new SystemEquetions(_interType);
            }
        }

        public static void init(stWindow stWin)
        {
            _stWindow = stWin;
        }

        public void initDefault()
        {
            longrange = defaultLongRange;
            interType = Routine.defInterType;
        }

        public static double defaultLongRange = 1.0;
        public double longrange;

        public double boxL;

        public int deg = 3;

        public MechSys(Atom[] atoms)
        {
            initDefault();
            foreach (var atom in atoms)
                if (atom != null)
                {
                    Add(atom);
                    deg = atom.deg;
                }
        }

        public MechSys()
        {
            initDefault();
            //Surprisingly there is nothing to be here...  
        }

        public void recognizeAtomInteractions()
        {
            foreach (Atom a_base in this)
            {
                a_base.AtomsToInteractWith = new List<object>();
                foreach (Atom a_inter in this)
                {
                    if (a_base == a_inter)
                        continue;
                    if (a_base.distance(a_inter) <= longrange)
                        a_base.AtomsToInteractWith.Add((object)a_inter);
                }
            }
        }

        public void addAtoms(Atom[] atoms)
        {
            foreach (var atom in atoms)
                if (atom != null)
                {
                    Add(atom);
                    deg = atom.deg;
                }
            recognizeAtomInteractions();
            int i = 0;
            foreach (Atom a in this)
            {
                atoms[i] = a;
                a.potE = _eq.Edw;
                i++;
            }
            _eq.atoms = atoms;
        }

        public double Mass()
        {
            double Mass = 0.0;
            foreach (var atom in this)
                Mass += atom.mass;
            return Mass;
        }

        public double[] CentrMass()
        {
            double[] CM = new double[deg];
            double M = Mass();

            for (int i = 0; i < deg; i++)
            {
                foreach (var atom in this)
                {
                    CM[i] += atom.mass * atom.pos[i];
                }
                CM[i] /= M;
            }

            return CM;
        }

        public double[] VelCenterMass()
        {
            double[] CMV = new double[deg];
            double M = Mass();

            for (int i = 0; i < deg; i++)
            {
                foreach (var atom in this)
                {
                    CMV[i] += atom.mass * atom.vel[i];
                }
                CMV[i] /= M;
            }

            return CMV;
        }

        public double[] Momentum()
        {
            double[] Mom = new double[deg];
            if (deg == 3)
            {
                foreach (var atom in this)
                {
                    Mom[0] += atom.pos[1] * atom.vel[2] - atom.pos[2] * atom.vel[1];
                    Mom[1] += atom.pos[2] * atom.vel[0] - atom.pos[0] * atom.vel[2];
                    Mom[2] += atom.pos[0] * atom.vel[1] - atom.pos[1] * atom.vel[0];
                }
            }
            return Mom;
        }

        public void MoveTo(double[][] vects)
        {
            int i = 0;
            foreach (var atom in this)
            {
                atom.MoveTo(vects[i]);
                i++;
            }
        }

        public double KineticEnergy(double[] vel)
        {

            double KE = 0.0;
            double KE0 = 0.0;
            foreach (var atom in this)
            {
                KE0 = 0.0;
                for (int i = 0; i < atom.deg; i++)
                {
                    KE0 += (atom.vel[i] - vel[i]) * (atom.vel[i] - vel[i]);
                }
                KE += 0.5 * KE0 * atom.mass;
            }
            return KE;
        }

        public double PotentialEnergy(Atom atom0)
        {
            double PE = 0.0;
            foreach (var atom in this)
                PE += atom0.Potential(atom0.Distance(atom.pos), atom.sort);
            return PE;
        }

        public double SelfPotentialEnergy()
        {
            double PE = 0.0;
            foreach (Atom atom1 in this)
                foreach (Atom atom2 in this)
                {
                    if (atom1.Distance(atom2.pos) != 0) PE += atom1.Potential(atom1.Distance(atom2.pos), atom2.sort);
                }
            return PE;
        }

        public void DoOneTimeStep()
        {
            double[] X = _eq.transformAtomsToSysVect();
            X = OdeSolver.OneStep(_eq.ODESys, _currT, _dt, X, null, IntegrationMethod.RungeKutta4);
            _eq.transformSysVectToAtoms(X);
            for (int i = 0; i < _eq.atoms.Length; i++)
                this[i] = _eq.atoms[i];
            _currT += _dt;
        }
    }

    /// <summary>
    /// The main class to interact with external programs
    /// </summary>
    public class Routine
    {
        public static sysTypes defInterType = sysTypes.Morze; 

        private StatusWindow _stWindow;
        public MechSys mSys;

        public Routine()
        {
            _stWindow = new StatusWindow();
            Parser.init(delegate() { return stWindow; });
            Atom.init(delegate() { return stWindow; });
            MechSys.init(delegate() { return stWindow; });
            _stWindow.Show();
            _stWindow.Disposed += new EventHandler(_stWindow_Disposed);
        }        

        public void parseInit()
        {
            mSys = Parser.readICfile(Parser.defaultFile);
        }

        #region status window handling
        void _stWindow_Disposed(object sender, EventArgs e)
        {
            stWindow = null;
        }

        public StatusWindow stWindow
        {
            get
            {
                if (_stWindow == null)
                {
                    _stWindow = new StatusWindow();
                    _stWindow.Disposed += new EventHandler(_stWindow_Disposed);
                    _stWindow.Show();
                }
                return _stWindow;
            }

            set
            {
                if (value == null)
                {
                    if (!_stWindow.IsDisposed && !_stWindow.Disposing)
                        _stWindow.Dispose();
                    _stWindow = null;
                }
            }
        }
        #endregion

    }
}
